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SEISMIC INVERSE SCATTERING IN THE ‘WAVE-EQUATION’ APPROACH 


CHRISTIAAN C. STOLK AND MAARTEN V. DE HOOP 

Abstract. Seismic data are commonly modeled by a high-frequency single scattering 
approximation. This amounts to a linearization in the medium coefficient about a smooth 
background. The discontinuities are contained in the medium perturbation. The wave 
solutions in the background medium admit a geometrical optics representation. Here we 
describe the wave propagation in the background medium by a one-way wave equation. 
Based on this we derive the double-square-root equation, which is a first order pseudodif¬ 
ferential equation, that describes the continuation of seismic data in depth. We consider 
the modeling operator, its adjoint and reconstruction based on this equation. If the rays in 
the background that are associated with the reflections due to the perturbation are nowhere 
horizontal, the singular part of the data is described by the solution to an inhomogeneous 
double-square-root equation. We derive a microlocal reconstruction equation. The main 
result is a characterization of the angle transform that generates the common image point 
gathers, and a proof that this transform contains no artifacts. Finally, pseudodifferential 
annihilators based on the double-square-root equation are constructed. The double-square- 
root equation approach is used in seismic data processing. 


1. Introduction 

In reflection seismology one places point sources and point receivers on the earth’s sur¬ 
face. The source generates acoustic waves in the subsurface, that are reflected where the 
mediumproperties vary discontinuously. The recorded reflections that can be observed in 
the data are used to reconstruct these discontinuities. 

The data are commonly modeled by a high-frequency single scattering approximation. 
This amounts to a linearization in the medium coefficient about a smooth background. The 
discontinuities are contained in the medium perturbation [|l]]. Thus a linear operator, the 
modeling operator, depending on the background, that maps the perturbation to the data 
is obtained. Both the smooth background and the perturbation are in general unknown 
and have to be reconstructed jointly. In this paper we analyze this reconstruction in the 
wave-equation approach. 

The reconstruction of the perturbation given the background is essentially done by ap¬ 
plying the adjoint of mentioned linear map (seismic imaging). The solutions in the back¬ 
ground medium admit a geometrical optics representation. Thus the modeling operator is 
Fourier integral operator [jT5p (for a general reference of Fourier integral operators see [§]). 
If the composition of adjoint and modeling operator, the normal operator, is pseudodiffer¬ 
ential, then the position of the singularities of the perturbation are recovered by applying 
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the adjoint to the data, and a microlocal reconstruction can be carried out. Under vari¬ 
ous assumptions on the background, concerning the presence of caustics and the geometry 
of the rays, results concerning the normal operator have been obtained 0,0.000] 
(section |2j.) 

In the Kirchhoff approach an approximation to the adjoint operator is constructed that 
depends on the background only through travel times and complex amplitudes that appear 
in the geometrical optics approximation. In the so called wave-equation approach 0, ^j, [T3|] 
the data are downward continued, leading to data from fictitious experiments below the 
surface at various depths. To form the adjoint a restriction is applied to the downward 
continued data (imaging condition). Using a one-way wave equation (section |4|) for the 
propagation of waves in the background, the downward continuation is described by the so 
called double-square-root equation. 

The data is formally redundant. There exist invertible Fourier integral operators that 
generate a set of reconstructions [ |I~8| ] from which the common image point gathers are 
obtained. Under Beylkin’s conditions 0], in particular in the absence of caustics, this can 
be done by using different subsets of the data. In the presence of caustics the mentioned 
set can be parameterized by angle between in- and out-going rays at the image point. If the 
background medium is correct the reconstructions in the set should be the same. This is a 
criterion that is used in the reconstruction of the background (migration velocity analysis). 
The redundancy leads to the existence of pseudodifferential operators that annihilate the 
singular part of the data [ fT8| J. 

In the Kirchhoff approach common image point gathers parameterized by angle can be 
generated by a generalized Radon transform. In the presence of caustics artifacts were 
observed in numerical examples in 0], By microlocal analysis of the Kirchhoff approach 
the presence of artifacts was shown in [jT7i]. In section [J] an approach to suppress these 
artifacts is discussed. 

In this paper, we consider the modeling, adjoint and reconstruction based on the double- 
square-root equation approach. The double-square-root equation is a pseudodifferential 
equation. If the rays in the background that are associated with the reflections due to the 
perturbation are nowhere horizontal, the singular part of the data is described by the solu¬ 
tion to an inhomogeneous double-square-root equation (section [5j). We derive a microlocal 
reconstruction equation in Proposition |6. I| . The main result, Theorem [7j] and Proposi¬ 
tion |7.2| . is a characterization of the angle transform that generates the common image 
point gathers, and a proof that this transform contains no artifacts. Annihilators based on 
the double-square-root equation are constructed in Corollary [8T| . 


2. High-frequency Born modeling and imaging 

We consider the scalar wave equation for acoustic waves in a constant density medium 
in K n . In preparation of the later analysis, we distinguish the vertical coordinate z G M 


pn— 1 


from the horizontal coordinates x G 
scalar acoustic wave equation is given by 


and write 


In these coordinates the 


Pu = f , P = c(x, z ) 


-2 


d_‘ 

m 


71—1 

£ 

3 = 1 


Dl+Dl 
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where D x = —i -§- c ,D z = —i Jr- The equation is considered in a time interval ]0, T[. 

If c G C°° the solution operator of ([[]) propagates singularities along bicharacteristics. 
These are the solutions of a Hamilton system with Hamiltonian given by the principal 
symbol of P 

P{x, Z, f, c, T) = -c(x, z)~ 2 r 2 + ||£|| 2 + C 2 . 


The Hamilton system is given by 

d(x,z,t) = dP d(£,C,r) _ dP 

d\ dX d(x,z,t)' 

Its solutions will be parameterized by initial position (x 0 , To), take-off direction a G .S' 7 '" 1 
and frequency r, 

X = x(xq, Zq, a, T, t) 


and similarly for z,t, r is invariant along the Hamilton flow. Here the evolution 
parameter is the time t. In section [f] we will change the evolution parameter to z, and use 
a similar notation to denote the bicharacteristics parameterized by z. 

By Duhamel’s principle, a causal solution operator for the inhomogeneous equation ([[]) 
is given by 


(3) 


u(x, z, t ) 



G{x, z, t - t 0 , x 0 , Zo)f(x 0 , Zq, t 0 ) d^od^dfo, 


Jo J 

where G is a Fourier integral operator with canonical relation that is essentially a union of 
bicharacteristics. Its kernel can be written as a sum of contributions 


(4) G{x,z,t,x 0 ,Zo) = V' / a^\x, z, t, Xq, Zo, 6) exp[i0^)(a:, z, Xq, Zo, t, 6)] dd, 


where the 0 d) are non-degenerate phase functions and the a® suitable symbols, see [§, 
chapter 5]. 

We adopt the linearized scattering approximation. The linearization is in the coefficient c 
around a smooth background Co, c = Co + Sc. The perturbation Sc may contain singularities. 
We assume that its support is contained in z > 0. The perturbation in G at the acquisition 
surface z = 0 is given by 


(5) SG(r, 0, t, s, 0) 



G(r, 0 , t 


to, To, Zo) 2 C 0 3 (x 0 , Zq)Sc(x 0, Zo) 

d 2 Q G(x 0 , To, to, s, 0) dtod^odTo, 


where both r, s G K n 1 . The singular part of SG is obtained by substituting (j?) into (5|). 
This defines the data modeling map 


F = F[c 0 ] : Sc^UoSG, 


where TZq is the restriction defined by 

TZq : I?p(M 2ri+1 ) -> V\Y) , u(x, z,t,xo, Zo) ^ (TZ 0 u)(x 0 , x,t) — u(x,0,t,x o ,0) 

with acquisition manifold Y a bounded open subset of R 2n-2 x M + that contains the 
range of values of (s, r, t). The restriction is defined on distributions V' v with wavefront 
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sets contained in an open conic set T such that T does not contain points of the form 
(x, 0, t , x 0 , 0, 0, C, 0, Co, 0). Since r ^ Oon WF(5G) the restriction is well defined. Since 
Y is bounded and the waves propagate with finite speed we may assume that Sc is supported 
in a bounded open subset X of M™” 1 x M + . We assume that X n {z = 0} = 0. 

Assumption 1. There are no rays from (s, 0) to (r, 0) with travel time t such that (s. r , t) G 
Y. For all ray pairs connecting (■ r, 0) via some (. x, z ) G X to (s. 0) with total time t such 
that (s, r, t) G Y, the rays intersect the plane z = 0 transversally at r and s. 


Theorem 2.1. [ [1~5[ |T0| l With Assumption 0 the map F is a Fourier integral operator V (X) 
—> T>'{Y) of order (n — l)/4 with canonical relation 


( 6 ) 

{ (x(x, z , p, T, t B ), x{x, z , a , T, t z ),t B + t z , £(x, z, f3, r, t B ), £(x, z, a, r, t z ), r; x, z, f, C) I 
t B , tr > 0, z(x, z, p, T, t B ) = z(x, z, a, T, t r ) = 0, (C, C) = -TCo(x, x) _1 (a + P), 

(x, z,a, P,t) G subset of X x (S n ~ 1 ) 2 x M\0}. 

Assumption |I] is microlocal. One can identify the set of points (s, r, t, a, p, r) G T*Y\ 0 
where this assumption is violated. If the symbol = ip(s,r,t,cr, p,r) vanishes on a 
neighborhood of this set, then the composition tpF of the pseudodifferential cutoff w = 
-0(s, r, t, D s , D r , D t ) with F is a Fourier integral operator as in the theorem. 

We assume that is as before and in addition vanishes outside Y. To image the singular¬ 
ities of Sc from the data we consider the adjoint F*i/j, which is a Fourier integral operator 
also. 


Assumption 2. ([7]] The projection of the canonical relation (|6|j on T*Y\ 0 is an embed¬ 
ding. 


Since @ is a canonical relation that projects submersively on the subsurface variables 
(x, z, C, C), the projection of (Q) on T*Y\ 0 is immersive (^. 25.3.6]. Therefore only the 
injectivity in the assumption needs to be verified [|T0|] . 

The following theorem describes the reconstruction of Sc modulo a pseudodifferential 
operator with principal symbol that is nonzero at (x,z,£, C) whenever there is a point 
(s, r, t, a, p, t; x, z,^,Q i n the canonical relation (g) with (s, r, t, a, p, r) in the support of 
(i.e. whenever there is illumination or insonification). 


Theorem 2.2. With Assumption |2] the operator F*v:F is pseudodifferential of order n — 1. 


3. Generalized Radon transform in scattering angle 


Consider the projection of the canonical relation <Q) on the (x, z, s , r, r) variables. Where 
this projection is locally diffeomorphic, the canonical relation (||) can be described by a 
phase function of the form 

r(T^ m \x, z, s, r ) — t) 


where is the value of the time variable in (|6). There is a set that de¬ 

scribes the canonical relation except for a neighborhood of the subset of the canonical 
relation where mentioned projection is degenerate. Each T (m) is defined on a subset D irn) 
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of M'^'7 2 r y We define F^ to be a contribution to F with phase function given by T( m \ 
such that on a subset of the canonical relation where the projection is nondegenerate F is 
given microlocally by M F^ n \ 

We can use (x, z, £, £) £ T*(M n_1 xR + )\0 as local coordinates on the canonical relation 
(6j). In addition, we need to parameterize the subsets of the canonical relation given by 
(x, z,£,() — constant; we denote such parameters by e. The canonical relation (Q) was 
parameterized by ( x , z, a, 3. r). We relate (x, z. £, e) by a coordinate transformation to 
(x, z, a, (3,t): A suitable choice when a 7 ^ /3 is the scattering angles given by 


(7) e(x, z , a, 3) = I arccos(a • /?), 


—01 T 3 


2 sin(arccos(a -3)/ 2) 


e]o, 


n —2 


The migration dip v is defined as 

(8) v( a >3) = 11° t All € S ‘ 

|| a + 3 1 | 


On there is a map (x, z, a , 3) | —► (x, z, s, r). We define (x, z, s, r) as the 

composition of e with the inverse of this map. Likewise, we define Y™') = Y m)> (x, z, s, r). 

We define the generalized Radon transform in scattering angle or the Kirchhoff angle 
transform via a restriction in F* of the mapping e3 m> to a prescribed value e, i.e. the 
distribution kernel of each contribution F ^ is multiplied by S(e — e( m )(x, z, s,r)) . Its 
kernel is given by 


(9) L(x, z, e, r, s,t) = YY~ [n ~ 1] / A (m) (x, z, s, r, T ) e ^ (7n H^,s,r,t, £ ,r) dr(k; 


meM 


where is a symbol for the m-th contribution to F, supported on I ) {rn) , and 

z, e, s, r, t, e, t) = T(Y m \x, z, s, r) — t) + (e, e — e^(x, z, s, r)). 


Here, e is the cotangent vector corresponding to e. Let 3>l = D r , D t ) be a pseudo- 

differential cutoff such that ij)(cr, p, r) = 0 on a conic neighborhood of r = 0. Then [[17!] 
3lL is a Fourier integral operator with canonical relation 


(10) U meM {(a;, e (m) (x, z, s, r), £ (m) (x, z, s, r, r, e), G (m) (x, z, s, r, r, e), e; 

s, r, T (m) (x, z, s, r) , p (m) ( x,z,s , r, r, e), cr (m) (x, z, s, r, r, e), t) \ 

(x,z,s,r) E D {m \ e e M" -1 , r G M\0} 


where 

(11) £ (m) (x,z,s,r,r,e) = d x $ {m) = Td x T {m) (x,z,s,r) - (£,d x e^ m) (x,z,s,r)), 

and likewise expressions for Y m> ■ an d p irn> . 

Let d be the Bom modeled data. To reveal any artifacts generated by L, i.e. singularities 
in Ld at positions not corresponding to an element of WF(5c), we consider the composition 
L F. This composition is equal to the sum of a smooth e-family of pseudodifferential 
operators and, in general, a non-microlocal operator the wavefront set of which contains 
no elements with e = 0 [if7|, theorem 6 . 1 ]. 
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We discuss the practical considerations of the suppression of artifacts generated by L by 
modification of transform (|9[). The correct way to remove artifacts from Ld is by a Fourier 
cut-off in e/||(£, (')| about 0. For bandlimited data, this is done approximately by local 
averaging in e at every (x. z ). Such procedure has been applied to both synthetic and real 
data examples in [fl|]. The artifacts were not fully suppressed. The following modification 
of the kernel of L (Q) was applied, 


L x (x,z,e,r,s,t) = ^ (2tt) (n x) 

m^M 


A( m \x, z, s, r, r)x(x, z, v^ m \x, z, s, r )) 

x e i*("‘)(x, Z , e ,a,r,t,e,A drd£; 


where x{ x i z i u ) is a smooth cutoff function on M n x S n 1 . Observe that A rn ) is the 
direction of a T ^ . 

Remark 3.1. The transform L x restricts the wavefront set of operator L. If WF(<5c) is 
contained in {(x, z, \{a x {x , z), a z (x , z))) | ( x , z, A) G M” _1 xl + x M}, where ( a x (x , z ), 
a z (x, z)) is a smooth covector field on M ra_1 x M+, there is a x with a small support in v 
such that L x generates the true image. Artifacts with singular direction such 

that the direction of is outside the support of %, are suppressed. 

The main difficulty arises when the background medium is not (accurately) known. 
Without knowledge of the background medium there is no criterion to distinguish arti¬ 
facts from the true image. The approach following the double-square-root equation to be 
introduced below does not generate artifacts. This is of particular importance in the devel¬ 
opment of tomographic methods for reconstructing the background medium. 


4. The one-way wave or single-square-root equation 


In this section we discuss the problem of solving the wave equation by evolution in the 
vertical, z direction. This problem is in general not well posed, but microlocal solutions 
can be obtained. 

Consider the wave equation rewritten as a first-order system 


( 12 ) 


dz {|) {~A(x,z,D„D t ) o)(|) + (/)’ 


where A(x , z, r) = c 0 (x, z) 2 r 2 — ||£|| 2 . Microlocally, away from the zeroes of A(x, z, 
A, r), this system can be transformed into diagonal form modulo a smoothing operator 
m- There is a family of pseudodifferential operator matrices Q(z) = Q(x, z. D Xl D t ) 
such that, microlocally, 



satisfy the one-way wave or single-square-root (SSR) equations 
(13) (-§z±iB±(x,z, D x ,D t ))u± = f±, 
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see [pG|]. The principal symbol b of the B± are given by b(x, z, £, r) = -\/A{x, z, £, r) = 
r _2 ||£|| 2 . For (x, t, £, r) such that the symbol B± is real, the equation is 


co{x,z) 2 


of hyperbolic type, corresponding (microlocally) to propagating waves. To describe the 
associated bicharacteristics the Hamiltonian P can be used, as well as (-fb(x , z, £, r). The 
solution of a one-way wave equation describes the propagation of singularities along rays 
in intervals where ±^| > 0. The case where the symbol B± is imaginary corresponds to 
either evanescent waves, or waves that blow up like in a backward heat equation. 

We choose the normalization of Q(z), such that (|T3j) is selfadjoint microlocally where 
the symbol is real, 


u = Q* + u+ + Q*_u_, 

f± = ± |i Q±f, 

where Q± = Q±(z ) = Q±{x, z, D x , D t ) are ^-families of pseudodifferential operators 
with principal symbols t~ 1/2 ( co ^ z)2 - 2 ||£|| 2 ) 1/4 . 

At the zeroes of A(x, z, £, r) the operators B±(z), Q±{z ) are not yet defined. To this 
end, we regularize the problem by replacing A(x, z , £, r) in ( fT2| ) by 

(14) c 0 (x, z)~ 2 r 2 - ||£|| 2 - ir 2 0(x, z, f, r), 


where <j>(x, z, r) is positive, small, homogeneous of order 0 in (£, r) and is supported on 
a small neighborhood of the set of zeroes of A(x, z, £, r). With this modification the opera¬ 
tors B±, Q± are ^-families of pseudodifferential operators, defined on the entire cotangent 
space T*W^ xt Y 0. With this choice of sign for the regularizing imaginary term there is a 
well defined solution operator G_(z, z 0 ), z 0 > z, of the initial value problem for ?j_ given 
by ( [Q| ) with /_ = 0, see [|22| . theorem XI.2.1]. The adjoint G-(z, z 0 )* describes the prop¬ 
agation from z to zq of (]13|), regularized in accordance with ([14|), but with opposite sign 
of the imaginary part. The operator G- is a Fourier integral operator with complex phase 
o a chapter XXV], [ 22] , chapters X and XI]. By Duhamel’s principle a microlocal 
solution operator for the inhomogeneous equation is given by 


(15) «_(-,z)=/ G_(z,z 0 )f_(-,z 0 )dz 0 . 

J — OO 

Microlocally —^i Q*_(z)G_(z, z 0 )Q-(z 0 ) is equal to the Green’s function for singularities 
propagating along bicharacteristics (cf. ([2])) with — > e, for some e > 0 that depends on 

c 0 and the support of <f> in (|T4l). 

The operator G_ propagates singularities at (x 0 , ( 0 , r, z 0 ) along the bicharacteristics for 
z in an interval \Z(x {) . ; ^ 0 , r, z 0 ), z 0 ], which is the maximal interval such that the regularized 
symbol b is real valued. As a consequence the bicharacteristic in this interval is nowhere 
horizontal. From now we use 2 as the evolution parameter for bicharacteristics, and denote 
them as 


(16) (®(X 0 , Zq, f 0 , T, z), Z, t(x 0 , Zq, f 0 , T, z), £(x 0 , Zq, f 0 , T, z), 

b(x(x 0 , Zq, ^0, T, z), Z, £(x 0 , Zq, ^ 0 , T, z),t),t). 
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Remark 4.1. The symbols of the operators B± (z ), ()± (z) can be written as an infinite sum 
of elementary symbols, Y^Li v i{ x i z i T ) w i{ z i^i T ) sa Y> that is rapidly converging where 
they are smooth, i.e. away from the set of zeroes of A(x, z, £, r), in accordance with [|Xf 


equation (2.1.11)]. The elementary symbols correspond to multiplications either in hori¬ 
zontal position space (iy) or in horizontal wavenumber (Fourier) space (w,). This is ex¬ 
ploited in fast numerical solvers of (|T3|), see [0. 


5. Modeling, the double-square-root equation 


We show that the Bom modeling operator can be written, modulo smoothing terms, 
in terms of the solution operator to the double-square-root (DSR) equation (21) below. 
We assume that the rays that connect source and receiver to a scattering point in X have 
nowhere horizontal tangent directions. 

Assumption 3. (DSR assumption) If (x, z) G X and a, 6 G S n ~ t s , t T > 0 depending on 
(x, z, a, (3) are such that z(x, z, /3, r, t s ) = z(x, z, a, r, t T ) = 0 and 
(x(x, z, (3, t, t s ), x(x, z, a , r, t T ),t s + t r ) G Y (cf. (§)), then 

dz 

— (x,z,f3,T,t) < —e, t G [0 ,t s \, 
dz 

— ( x,z,a,r,t) < —e, t G [0,t r ], 

where e > 0 was introduced below 0- 

This assumption is stronger than Assumption [[[ This assumption is microlocal, and, 
given the background medium, a pseudodifferential cutoff can be applied to the data to 
remove microlocally the part of the data where Assumption |3] is violated. 

Under Assumption [T] and the assumption that dc = 0 on a neighborhood of l = 0, the 
singular part of the Bom data is unchanged when G in (^j) is replaced by 
—|i Q*-(x, z, D x , D t )G-Q-(x , z, D x , D t ). Define the operator H(z, zq), z < z 0 by 


(17) (H(z,zo))(s,r,t,s 0 ,r 0 ,to) = 


(G-(z, z 0 ))(s, t — t 0 


t', s 0 , 0) (G- (z, z 0 ))(r, t\ r 0 , 0) df ; . 


Here (G_(z, z 0 ))(r, t'. r 0 , 0) denotes the distribution kernel of G_(z, z 0 ), and similarly for 
H{z, ^o)- Define the maps X\ , X 2 by 


(18) Ji : D'(R n ) -h. ©'(M 2 ”" 1 ) : u(x, z) ^ S(r - s)u(z), 

(19) X 2 : D^R 2 "- 1 ) ^ D'(R 2n ) : u(r, s, z) ^ 5(t)u( r -± s -, z). 

The operators G, G_, Q_ all are of convolution type in the time variable. It follows that 
the singular part of the Bom approximated data (5) is given by 


max 


(20) FSc = / Q*_ a (0)Q*_ r (0)H(0, z)Q^ s (z)Q., r (z) \D 2 t { X 2 X x cf5c)(-, z) dz. 

Jo 

where Q_ >s (z) is short for Q-(s, z, D s , D t ) and similarly for Q^^ r (z). Here z max is the 
maximum depth illuminated from acquisition manifold Y given the background medium 


Co- 
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Define the inhomogeneous double-square-root (DSR) equation by 

(21) (£ - \B_{s,z,D s ,D t ) - i B_(r,z,D r ,D t ))u = g(s,r,t,z). 

It follows from the definition of G-(z, z 0 ) that the operator H(z, z 0 ) is a solution operator 
for the Cauchy initial value problem for the (regularized) DSR equation. The Born approx¬ 
imated data is modulo smooth functions equal to the solution of the DSR equation at z = 0 
with inhomogeneous term g(s, r, t, z) = \D\(T 2 T\ Cq 3 <5c) (s, r, t, z) . 

The operator H(z, zq) propagates singularities along the bicharacteristics for z in an 
interval ]Z min (s 0 , r 0 , cr 0 , po, r), zq\, which is the intersection of the one-way intervals asso¬ 
ciated with the source and receiver bicharacteristics, in the notation of ( l_6j). 

(22) (x(so,zo,a o ,T,z),x(r o ,z o ,po,r,z),t 0 + t(s 0 , z 0 , a 0 ,T, z) + t(r 0 , z 0 , p 0 ,r, z), z, 

£(s 0 , Zo, *0, T, Z ), £(r o, Zq, po, T, z), T, 

r(*(s 0 , Zq, CTq, T, z),x{tq, Zq, p Q , T, z), £(s 0 , Zq, CTq, T, z), £{t Q , Zq, po, T, z), T, z)), 

where 

(23) r(s, r, a, p, r, z) = b(s , z, a, r) + b(r, z, p, r). 

6. DOUBLE-SQUARE-ROOT RECONSTRUCTION 

Let N := F*ipF, where ij) — ip(s, r, t, D s , D r , D t ) is a suitable pseudodifferential cutoff 
as in section By Theorem [2~Tj , with Assumption N = N(x, z, D x , D z ) is a pseudodif¬ 
ferential operator. Therefore 

(24) N(x, z, D x , D z )Sc = F*i>(s, r , t D s , T> r , A)d, 

and <5c can be reconstructed microlocally where the principal symbol of N is non-zero. 
Assumption j3] implies the injectivity in Assumption [| and is hence stronger than Assump¬ 
tion ^jby the remark below Assumption [|. We present a reconstruction formula similar to 
( p4| ) based on the DSR Bom modeling ([20]). This leads, again, to reconstruction modulo a 
pseudodifferential operator for which an explicit expression is given. 

The adjoint operator H( 0, z)* propagates the data downward and backward in time. We 
consider H(0,z)*d as a function of ( s,r,t,z ). The operator H(0,z)* is microlocally a 
Fourier integral operator with real phase and canonical relation which follows from (|22j) 
with z 0 = 0, for z in an interval [0, Z max (s 0 , r Q , a 0 , p Q , r) [ 

(25) {(®(s 0 , 0, cr 0 , r, z), x(r 0 , 0, p 0 , r, z),t 0 + t(s 0 , 0, a 0 , r, z) + t(r 0 , 0, p 0 , r, z), z, 

€(So, 0, (To, T, z), £(r 0 , 0, Po, T, z), T , r(s 0 , To, (To, po, T, z)\ S Q , T 0 , t Q , (T 0 , po, t) \ 

(s 0 , r o, to, (T 0 , po, t) e T*M 2n_1 \0, z G [0, Z max (s 0 , r 0 , a 0 , p 0 , r)[}. 

The adjoint of the operator J 2 (cf. ([T9[)) is given by the restriction 72. 2 defined by 

g(r,s,t,z) ^ ( n 2 g)(r,s,z ) = g(r,s,0,z). 

If u = u(s, r, z) then Ku = ( Ku)(s, r, t) will be defined by 

Ku — I H(Q, z)I 2 u(-,z) dz. 

J — OO 
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Let ip be a pseudodifferential cutoff in the (s, r, t) variables supported in 

{(s 0 ,r 0 ,t 0 ,a 0 ,p 0 ,T) 1 t 0 £ 

[0, t(s 0 , 0, (To, T, -^max(So, To, <To, Po, t)) tffoi 0, Po, T, ^max('Sc), To, (To, Po, t)) [^ . 

Then ?pK is a Fourier integral operator with real phase. The adjoint K* of K follows from 
the equality 

77(0, z)l 2 u(-,z)dz)( St r tt ) = f (U 2 H(0,z)*'iJ;d,u(‘,z))( Str )dz 

J —OO 

and is given by 



K*f = n 2 H{ 0, z)*-0. 

Since < 0, t depends strictly monotone on z for each (s Q , r 0 , t 0 , <r 0 , p 0 , r). The DSR 
bicharacteristic with initial values (s 0 , r 0 , to - <To, Po, t) hence intersects the /; = 0 hyper¬ 
plane at most once, and the intersection is transversal. From this, it follows that the com¬ 
position K*ip is a Fourier integral operators with a locally invertible canonical relation. 

The canonical relation maps points (s 0 , r 0 , t 0 , <r 0 , p 0 , r) in a subset of the cotangent ac¬ 
quisition space T*M 2n_1 \0 to (r, s, z, <r, p, C) in a subset of T*M 2n_1 \0. This map converts 
time to depth. We define the symbol T (s, r, z, a, p, £) as the pull back of ip(s 0 , r 0 , t 0 , cr 0 , 
Po, r) under the inverse of this map. Starting from the Born modeling (GO), we find the 
following reconstruction proposition. 

Proposition 6.1. There are pseudodifferential operators $ = z, D :r . D z ) of order 
7i — l with principal symbol 

( 26 ) $(x,z,£,()= [ d>(x,x,z, ^-9,^ + 9,C)d9, 

J R"- 1 

and S(z) = S(s, r, t, D s , D r , D t , z) of order 0 with principal symbol 


l(s,r,t,o,p,r,z ) = 


ar 

57 


(S, T, (T, p, T, Zj 


= c 0 (s,z) (c 0 (s,z) 


' ~ 2 -T- 2 


IM| 2 ) 1/2 + c 0 (t,z) 2 (c 0 (t,z) 


-2 _ T — 2 | U || 2 ^- l /2 


such that 

(27) $(x,z,D x ,D a )(5c 

= 2c^^ 2 s(z)g*_ iS (z)- 1 gi ir (z)- 1 //(0,z)*g_ iS (0)- 1 g_, r (0)- 1 A" 2 V’c(, 

where d = Fc)c is the Born modeled data. 

Proof We calculate microlocally the principal symbol of K*K. The kernel of the operator 
77(0, z) has microlocally an oscillatory integral representation with a phase function asso¬ 
ciated with generating function, S = S(z, s, r, t, y 0I , p 0 j) say, where p 0 = (so, t 0 , to) and 
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r/ 0 is the corresponding cotangent vector, and {I, J} is a partition of {1,... , 2n — 1}, 

H(0,z)(y 0 ,s,r,t) = (2n)-( 2n - 1+ W 2 

x J A(z, s, r , t, yo, r 10 j)e^ s ^ z,s,r ’ t ’ yoI ’ VoJ ^~^ VoJ,yoJ ^ dr/ 0 j, 

The adjoint H( 0, z)* has amplitude A(z, s, r, t, y 0 , yoj) and phase — 5( 2 , s, r, £, y 0 /, r] 0 j ) + 
(r/oj, yoj) • Hence, the kernel of the composition 7/(0, 2 )*7/(0, ;j) has the oscillatory inte¬ 
gral representation 

(27r) -(2n-1) J A(z, s', r', t', y 0 , r} 0J )A(z, s, r, t, y 0 , y 0J ) 

X e i [- s ( z > s '> r ’> t '>Voi,‘noj)+S(z,s,r,t,yoi,Voj)] 

We expand the phase in a Taylor series around (s', r', t') = (s, r, t ) and identify the gradient 

dS 

(a(z, s, r, t, y 0I , y 0J ),p(z, s, r, t, y 0I , y 0 j), r(z, s, r, t, y 0I , rj 0J )). 
Applying a change of variables, (y 0 /, r] 0 j) (—>• (a, p, r), the phase takes the form 

(fo p, r),(s'-s,r'-r, f'-f)). 

Since 7/(0, 2 )*77(0, z) is a pseudodifferential operator with symbol 1, microlocally in the 
support of the cutoff ip, we conclude that the principal part a of the amplitude A satisfies 

The kernel of operator K has an oscillatory integral representation similar to the one of 

H{0,z), 

K(y 0 ,s,r,z ) = (2tt )-V*-i+\iW J A ( z , s, r, 0, y 0 , 

(we have applied J 2 at z). It follows that the composition K*K, carrying out an analysis 
similar to the one for H{ 0, z)*H{ 0, z), is a pseudodifferential operator, microlocally. Its 
amplitude has principal part 

P ) _1 d(a,p,r) 
d(y 0 i,poj) d(y 0 i,r) 0 j) 

For fixed (s, r, a, p, z) the map r h-> T(s, r, a, p, r, z) is invertible on a set given by \t\ 
sufficiently large. This map will be denoted by T _1 = r _1 (s, r, z, a, p, Q. It follows that 
the principal part of K K is microlocally given by 

fjY 

— (s, r, a, p, r _1 (s, r, z,a,p,Q, z ) 


| a(z, s, r,t,y 0 , Vo j) \ = 


9((t,p,t) 

d(y 0 i,Voj) 



-1 
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We finally consider the composition of operators TZ\^X\. Its kernel has an oscillatory 
integral representation, 


(2vr)- (2n " 1) 
(2t r 


JR 2 ™- 1 

— (2n—1) I 


V(x, x, z, a, p, C) JSUwO) dpdadC = 

[ y(x,x,z,#-e,# + 9,0 d^e^-^+^-^dedC, 


upon changing variables of integration, a = — 6, p = + 9. The domain of the 

6 integral is bounded depending on (£, (), since 'k is a cutoff in (s, r, t, a , p, r). Hence 
TZ \'kXi is a pseudodifferential operator of order n — 1 with principal symbol ([26j). 

Applying the above results to the expression ([20]) for the Bom modeled data leads to the 
statement of the proposition. □ 


Remark 6.2. Note that in (27) all the operators on the right hand side, except the down¬ 
ward continuation 77, act at depth z only. On the contrary, the operator $(x, z, D X ,D Z ) 
on the left hand side depends on the Hamiltonian flow associated with the background 
medium in the depth interval [0, z]. In the usual wave-equation imaging algorithms it is 
hence straightforward to include the pseudodifferential factors on the right hand side of 
(27). To account for the operator <k(x, z, D x . D z ) on the left hand side one requires an 
additional ray computation. 


Remark 6.3. Depending on the background medium, the reconstruction can also be done 
using data on a submanifold Y' of Y. Let TV be the restriction of a function on Y to 
Y\ so that the forward map for this case is given by 7 Z'F. In suitable local coordinates 
(y f , y") on Y such that y" = 0 defines Y', the adjoint X' of 7 Z' is given by the map 
(X'f)(y',y") = f(y')8(y"). Conditions such that F*X ,/ ijj , TZ'F is pseudodifferential are 
given in 0, where , 0 / is a suitable pseudodifferential cutoff. Reconstruction modulo a 
pseudodifferential operator is done in this case by first applying the map X' to the data, and 
then applying the previous procedure. Applying X' to the data simply means adding zeroes 
where there is no data in Y. 

7. The wave-equation angle transform 

We define the wave-equation angle transform /1 WE by the following integral of the 
downward continued data, 77(0, z)*d, 

(28) ( A WE d)(x, z,p) = I (77(0, z)*^d)(x — \,x + %,ph)x(x, z, h) dh, 

J i"- 1 

(cf. 0), where h > y(a;, z, h) is a compactly supported cutoff function the support of 
which contains h — 0. 

Theorem 7.1. Let Co be an upper bound for c 0 . Assume that 

(29) \p\ < Pmax < \Cff l . 

Then A WE is a Fourier integral operator such that A WE F is a smooth p-family ofpseu- 
dodifferential operators in (./;, z). Let C\ be an upper bound for If in addition the 
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function h i—> x(x, z, h) is supported in B( 0, R), where R depends on Cj and 'if then the 
canonical relation of Awe corresponds to an invertible map from a subset ofT*M. 2 ™~^ to 
a subset of T* Mj-i) 

that has nonempty intersection with the set 9 = 0 (where 9 is the 

p-covector). 


Proof The map d i—> 11(0, z)*v>d is a Fourier integral operator with canonical relation 
given in (j25|). 

The Schwarz kernel of the map H( 0, z)*f>d i—> A WE d equals 


5(x — 5(p(r — s) — t) S(z — z') z,r — s) 

= (2n)~ n ~ 1 J t )+r(p(r-s)-t)+c(z-z')) d^drdC- 

It is a Fourier integral operator with canonical relation that is contained in T* >< 

™£r,M)\° and s iven b y 


(30) c, (r - s)r; s,r,p(r - s),z , § +pr, § -pr,r, () \ 

(s, r, z, p, £, (, t) G (subset of )M 4 ” -1 }. 


This canonical relation can be parameterized by the coordinates of T*R 2 ™ rtz ^\0 except t. 


that is (s, r, z, o, p,r,Q. The projection of ( 0O| ) on T* 
by 


o2 n 


s,r,t,z ) 


\0 is a hypersurface defined 


(31) 


t = 


°~P 
2 r 



The canonical relation of the map d i—► ff(0, z)*f>d , considered as a function of (s, r, t, z), 
is parameterized by (s 0 , r 0 ,t 0 ,o 0 , p 0 , r, 2 ). The canonical relation ( [25j ) is time translation 
invariant and the line in T*Mf ( f rt AO parameterized by t 0 for fixed (s 0 , r 0 , o 0 , Po, t) inter¬ 
sects the hypersurface ( pT| ) trans vers ally. It follows that the composition of the canonical 
relations ( |25| ) and ( ]3F| ) is transversal. The composition is parameterized by (s 0 , r o, cr 0 , p 0 , 
r, z). It follows that ylwF. is a Fourier integral operator. 

The composition if (0, z)*ipF, that maps 5c = 5c(x, z') to the downward continued data 
as a function of (s, r, t, z), is a Fourier integral operator with canonical relation 


(32) {(a;(a:, z! , o' , r, z), x(x, z', p' , r, z), t(x, z! , o' , r, z) + t(x, z', p' , r, z), z, 

z\ o', T, z), £(x, z', p, t, z),t,(;x, z', o' + p', (') I 

(. x,z',o',p'X',z ) G a subset of R 3n , r such that (' = T(x, x, o', p', r, z'), 

C = T(x(x, z', o', T, z), x(x, z', f>, T, z),£(x, z', o', T, z), £(x, z', p', T, z), T, z)}. 

The propagation of singularities upward by F and downward by ii(0, z)* is along the same 
DSR bicharacteristics. 

We show that the composition A WE F is a Fourier integral with canonical relation con¬ 
tained in 


(33) {(ar,2,p,f,C,0;a;,2,f,C) | (x,z,^,C) e T *^ X , Z )\Q,P e] - p max ,Pmax[}- 
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From that it follows that A^ E F is a p-family of pseudodifferential operators. The projec¬ 
tion of ( f32[ ) on T*Wiff r t ,. intersects the hypersurface ( pT| ) at z' = z, since then r — s = 0 
and t — 0, leading to elements in (03|). Since singularities propagate with speed less 
than C 0 , if (u s , v r , v t: . v z , v a , v p , 0, vfi is a tangent vector to the DSR bicharacteristic, then 
< 2C 0 . Therefore, by (|29l), the composition of ( |32| ) with ( ]30| ) is transversal and 
contains no elements outside 

The projection of the canonical relation of A WE on the second component T*Mf n ~\\ 0 is 
invertible if each DSR bicharacteristic with initial values (s 0 , r 0 , t 0 , a 0 , p 0 , r), parameter¬ 
ized by z, intersects the hypersurface d3l| ) at most once and trans vers ally. Let p(z) denote 
2^ along a certain DSR bicharacteristic and let t(z) denote the time and h(z) denote the 
value of r — s. The elements of the canonical relation of A WE correspond to solutions of 
t(z) — (p(z ), h(z)) = 0. To estimate the derivative of the left hand side we observe that 


dt dh 

Tz ~ {p(z) ’ a; (z)) K - £ " 

for some e 0 > 0 depending on the cutoff w (or <p in (|T4|)) and on the value p max . Since 

d£ db t 8cq 2 

dz dx ^Cq 2 - r -2 ||£|| 2 dx ’ 


there is a bound on ^ in terms of Cj and on u> (to bound the square root from below). It 
follows that for some 6\ < eo 


dp 

[ dz 


(z),h(z)) 


<e 1 <e 0 


if ||/r|| < C 2 Ci 1 for some constant C 2 depending on -0 and p, n:ix . This implies that the 
function z >—> t(z) — (p(z),h(z )) is monotone. Hence the projection of the canonical 
relation of A WE on 0 is invertible. It follows from the above reasoning that the 

projection on is invertible as well. This establishes the last statement of the 

theorem. □ 


To conclude this section we determine, at the principal symbol level, the modification of 
( |28l ) that leads to microlocal reconstruction. 

Proposition 7.2. Define /1 WE by 


(A WE d)(x,z,p) = / x{x,z,h)2c 0 (x,z) 3 

J R "- 1 

x (E(z)Q*_ s (z)~ 1 Q*_ r (zy 1 H(0, z)*Q_ tS (0y 1 Q_^(0y 1 Dfi 2 'fid) (x - 1, x+ 1, ph ) dh. 

Suppose that x(x, z, 0) = 1 and h i—> x(x, z, h ) is supported in B( 0, R ) (cf. Theorem m 
then A we is an invertible Fourier integral operator. Let the symbol T W e = 'ITveOt- z , p, 

C, d) (where 6 is the p-covector) be given by the pull back of w under the map from 
to T*M. 2 ff ^\0 induced by the canonical relation of A WE . The composition 

y/ l\VE F is a p-family of pseudodifferential operators with principal symbol 'I ; we(t, z,p, 

C,o). 
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Proof. It is sufficient to show that the map 


(34) 


5c 


i(z)H(0,z)* / H(0, z')IfL\5cdz' ) (x — x + |,p/i) dh, 


microlocally has principal symbol equal to 1. In this proof we will omit the cutoff functions 
that are part of the symbols; the calculations will be valid microlocally on the support of a 
cutoff. 

Using an oscillatory integral representation of H similar to the one in the proof of Propo¬ 
sition |5. ![ , we find that the principal contribution to the kernel of this map, as a function of 
(x, z, p; x', z 1 ), can be written as 


( 2*0 


dS 


( 2 n i) / s(x — Ah,x + \h,ph, ——, 


ds 


dS_ _dS_ . 

dr ’ dt ’ 


x A(z, x - \h, x + \h, ph, y 0 , t]oj)A(z', x', x', 0, y 0 , rj 0J ) 

X e i [- s ( z ’ x ~^ h ’ x +^ h ’P h ’yoi> r loj)+S(z',x’,x’,0,y O i,voj)] ^ y OI dl]oj dh 

We expand the phase in a Taylor series around (x', z', h) = (x. z, 0) and identify the gradi¬ 
ent at (x, z, 0), 


dS 

dx 

dS_ 

dz 

dS_ 

~dh 


(z, x, x, 0, y 0I , t)qj) = a{z, x, x, 0, y 0I , <q 0J ) + p(z, x, x, 0, y 0I , rj 0J ), 
(z, x, x, 0, y QI , r/ 0J ) = ((z, x, x, 0, y 0I , tj 0J ), 

(z, x, x, 0, y 0I , rjoj) 


= -\<t(z, x, x, 0, y 0I , p 0 j) + | p{z, x, x, 0, y 0I , r] 0J ) + pr(z, x, x, 0, y 0I , r] 0J ). 
Applying a change of variables, (y 0I , p 0 j) i—»■ (cr, p, (), the phase takes the form 

(cr + p,x - x') + C (z - z') + {\{p - a),h). 

The amplitude factor EAA becomes equal to one by the calculations in the proof of Propo¬ 
sition |6~i] Upon changing integration variables, a = — 9, p — + 9, the oscillatory 

integral takes the leading-order form (microlocally) 

( 27 r)- (2n - 1) J e K(£,*-*')+< {z-z')+(o,h)) d£d(d0dh. 


It follows by integrating the 6, h variables that ( |34j ) indeed has principal symbol equal to 
microlocally. 


1 

□ 


8. Annihilators 


It was observed in [[18]] that there are pseudodifferential operators that annihilate the 
data, due to the fact that the inverse problem is formally overdetermined. On transformed 


data A WE d, these are given by i = 1,... 


,n 


1 hence the annihilators are given by 


(-Awe) 


-i 


d_ 

dpi 


A 


WE, 





16 


CHRISTIAAN C. STOLK AND MAARTEN V. DE HOOP 


Where (^we) 1 is a regularized inverse for A WE , that is a microlocal inverse for a subset 
of T where A WE is invertible. 

We define the operator K* that maps d to a function (s, r, z) by 

(K*d)(s,r,z) = 2c 0 ((s + r)/2,z) 3 

x (n 2 E(z)Q*_ >s (z)- 1 Q*_ jr (z)~ 1 H(0, zyQ.^O^Q-A^D^d) (s, r, z). 

Using (]27|), we observe that the operator K* acting on d = FSc yields 

K*f>d = ^(s, r, z, D s , D r , D z )li5c 

(for the definition of T see remark above Proposition |67T| ). Applying the operator A4 given 
by the multiplication by r — s to this equation yields 

M.K*f>d = [Ad, v I ; (s', r, z, D s , D r , D z )] Ti Sc, 

i.e. there is only a lower order contribution, since (r — s)S(r — s) = 0 and hence A41 t = 0 
(cf. ©). 

For the subset in T*R^ r ^ where the operator fy(s,r, z, D s , D r , D z ) is microlocally 
elliptic, the operator 

M : = M - [M,^(s,r, z, D s , D r , D z )]^(s,r, z, D s , D r , D z )~ x 

is an annihilator of K*wd to all orders. 

Corollary 8.1. A pseudodifferential annihilator of the data is given by 

W = KMK*. 

Note that W = W [co] depends on the background medium. The semi-norm ||lU[co]d|| 
can be viewed as the wave equation analog of the differential semblance functional [ 19]. 
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